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ABSTRACT 

Spectral methods are well suited for solving hydrodynamic problems in which the self-gravity of the flow 
needs to be considered. Because Poisson’s equation is linear, the numerical solution for the gravitational poten¬ 
tial for each individual mode of the density can be pre-computed, thus reducing substantially the computational 
cost of the method. In this second paper, we describe two different approaches to computing the gravitational 
field of a two-dimensional flow with pseudo-spectral methods. Eor situations in which the density profile is in¬ 
dependent of the third coordinate (i.e., an infinite cylinder), we use a standard Poisson solver in spectral space. 
On the other hand, for situations in which the density profile is a delta function along the third coordinate (i.e., 
an infinitesimally thin disk), or any other function known a priori, we perform a direct integration of Poisson’s 
equation using a Green’s functions approach. We devise a number of test problems to verify the implemen¬ 
tations of these two methods. Einally, we use our method to study the stability of polytropic, self-gravitating 
disks. We find that, when the polytropic index T is < 4/3, Toomre’s criterion correctly describes the stability 
of the disk. However, when T > 4/3 and for large values of the polytropic constant K, the numerical solutions 
are always stable, even when the linear criterion predicts the contrary. We show that, in the latter case, the 
minimum wavelength of the unstable modes is larger than the extent of the unstable region and hence the local 
linear analysis is inapplicable. 

Subject headings: accretion disks — black hole physics — hydrodynamics 


1. INTRODUCTION 

In the standard model of accretion disks, turbulent viscosity 
plays an important role in bringing materi al inward and trans - 
porting angular momentum outward (see iFrank et alJl2002h . 
At the same time, viscous dissipation converts gravitational 
potential energy to thermal energy and heats up the disk, 
which then radiates away this energy as thermal emission. In 
most applications of accretion disks around central objects, as 
in, e.g.. X-ray binaries, the self-gravity of the flow is negligi¬ 
ble compared to that of the central object. However, there are 
disk-like systems, such as active galactic nuclei as well as pro- 
tostellar and protoplanetary disks, where the effects of self¬ 
gravity c hange not only the properties of angular momentum 
transp ort jBossll998HBalbus & Papaloizouli99^lMeifa et dJ 
20051). but also the energy b alance equation jBertinlll99'f 
Bertin & Lodatolll999t 1200 Ih . which affect the global struc¬ 
ture of the disk. 

Besides the angular momentum transport problem, self- 
gravitating disks are also very important in studying star and 
planet formation. Indeed, gravitational instabilities in proto¬ 
planetary disks have been proposed as viable planet forma¬ 
tion mechanisms. Although there has been a large amount of 
work done o n gravitational instabilities in these sys tem (see, 
for example. iPickett et alJl200^ iMeifa et alJl200-5[ and ref¬ 
erences therein), it still remains an open question whether 
the fragmentation by gravitational instabilities can produce 
bound planetary objects, or if reservoirs of small solid cores 
are needed for the accretion mechanism to lead to rapid planet 
formation. 

In the first paper in this series (Chan, Psaltis, & Ozel 2005), 
we presented a pseudo-spectral method for solving the equa¬ 
tions that describe the evolution of two dimensional, viscous 
hydrodynamic flows. There, we addressed issues related to 
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the implementation of boundary conditions, spectral Altering, 
and time-stepping in spectral methods, and verified our algo¬ 
rithm using a suite of hydrodynamic test problems. In this 
second paper of the series, we present our implementation of 
a Poisson solver that allows us to take into account the effects 
of self-gravity of the flow. 

Spectral methods are particularly suitable for incorporating 
the effects of self-gravity. Because Poisson’s equation is lin¬ 
ear, we can pre-compute the numerical solution for the gravi¬ 
tational potential for each individual mode of the density, thus 
reducing substantially the computational cost of the method. 
In fact, Fourier methods have been used extensively in ana¬ 
lytical studies of gravitati onal potentials in systems w ith peri¬ 
odic boundary conditions dBinnev & Tremaineil987h . On the 
other hand, in numerical studies of self-gravitating disks, the 
strengths of spectral methods have only been partially incor¬ 
porated. 

Hybrid hydrodynamic algorithms with self-gravity have 
been developed in such a way that modified spectral methods 
are used only for solving Poissson’s equation for the gravita¬ 
tional held, whereas the hydrodynamic parts are still treated 
with fi nite difference schemes. For example, iBoss & Mvhilll 
dl992h describe a spherical harmonic decomposition method 
and a second-order scheme i n the radial direction t o solve 
Poisson’s equation, whereas iMvhill & BossI use a 

modified Fourier method to And the gravitational potential of 
an isolated distribution of sources. Both algorithms employ 
explicit second-order finit e difference methods to ad vance the 
hydrodynamic equations. IPickett et alJ in^l2(i(i(il) describe 
an implementation of an algorithm that uses Fourier decom¬ 
position in the azimuthal direction of cylindrical coordinate 
to solve Poisson’s equation together with a von Neumann & 
Richtmeyer AV scheme for the hydrodynamics. Some other 
examples can be found in Grandclement et al. (2001), Brod¬ 
erick & Rathore (2004), and Dimmelmeier et al. (2005). The 
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main advantage of using hybrid methods is that one can em¬ 
ploy currently available hydrodynamic algorithms based on 
finite difference schemes. However, a hybrid method does 
not exploit the high order of the spectral algorithm, because 
the hydrodynamic difference schemes typically have a much 
lower order compared to that of the Poisson solver. Con¬ 
trary to these efforts, our algorithm uses a spectral decom¬ 
position method for solving both the hydrodynamic and Pois¬ 
son’s equation, providing a consistent treatment of the whole 
problem. To our knowledge, this is the first time that spectral 
methods have been used in studying astrophysical disks with 
self gravity. 

By construction, there is an ambiguity in the definition of 
the gravitational field in two-dimensional problems. We can 
assume either that the density profile is independent of the 
third coordinate (i.e., an infinite cylinder) or that it is a delta 
(or any other predetermined) function along the third coordi¬ 
nate (i.e., an infinitesimally thin disk); the resulting gravita¬ 
tional field on the two-dimensional domain of solution will 
be different in the two cases. For example, the gravitational 
potential of a “point source” on the two-dimensional domain 
of solution will be proportional to log(r) in the first case and 
to 1 jr in the second, where r is the distance from the source. 
In order to consider both geometries, here we describe two 
different approaches to computing the gravitational field of 
a two-dimensional hydrodynamic flow with pseudo-spectral 
methods. When the flow has the geometry of an infinite cylin¬ 
der, we use a standard two-dimensional pseudo-spectral Pois¬ 
son solver, which has been proven to be numerically stable 
and accurate. When the flow has the geometry of an infinitesi¬ 
mally thin disk, we perform a direct integration of the Green’s 
function for the gravitational potential, following the work of 
Cohl&Tohline(1999). 

In the following section, we present our assumptions and 
equations. In ^ we discuss the details of our numerical 
methods, include both the standard Poisson solver and the 
Green’s function integrator. Next, we present a series of tests 
in ^ to verify our algorithm. Finally, we apply our method 
to a numerical study of Toomre’s stability criterion of self- 
gravitating disks in §5. 

2. EQUATIONS AND ASSUMPTIONS 

We consider two-dimensional, viscous, compressible flows. 
In this second paper, we include self-gravity and continue to 
neglect the magnetic fields of the flows. The hydrodynamic 
equations, therefore, contain the continuity equation 

BT, 

—+ V.(Sv)=0, (1) 

the Navier-Stokes equation 
By 

E —-i-S(v V)v = -VP-i-VT-i-Eg, (2) 

and the energy equation 
dE 

-^-l-V • {Ey) =-PV •v-l-$-V -q- V-F-2F’j. (3) 

We denote by E the height-integrated density, by v the fluid 
velocity, and by E the thermal energy. In the Navier-Stokes 
equation, P is the height-integrated pressure, r is the viscosity 
tensor, and g is the gravitational acceleration. We use $ to 
denote the viscous dissipation rate, q to denote the heat flux 
vector, and F to denote the radiation flux on the r-cj) plane. 
The last term in the heat equation, 2F^, takes into account 


the radiation losses in the vertical direction. The analytical 
forms of the various physical q uantities in equations (□-(S 
are given in IChan et M ll200-5l) . except for the gravitational 
acceleration g, for which we need to introduce a new equation. 

In Newtonian gravity, the gravitational field g is conserva¬ 
tive, so we can define the gravitational potential T* by 

g = -VT-. (4) 

The gravitational potential associated with the (three- 
dimensional) mass density, p, is given by the volume integral 

over all space, where G is the gravitational constant. Rewrit¬ 
ing equation 0 in differential form, we obtain Poisson’s 
equation 

v2T-=47rG/9, (6) 

with 4* satisfying the boundary condition T'(f,oo) = 0 at all 
times. When simulating astrophysical flows, the computa¬ 
tional domain is usually finite. Based on its linearity, 
we can decompose the Poisson equation into two parts 

V^4'int = 47rGpint, (7) 

and 

V^4'ext = 47rGpext, (8) 

where pjnt denotes the mass density within the computational 
domain, which in our case is the flow density, and pext refers to 
external sources such as a central object and/or a companion 
star. The gravitational field is then given by 

8 — Sint + 8ext = ~^(4/int + 4'ext) ■ (9) 

In the astrophysical context of interest here, 4'ext is usually 
generated by a set of spherical objects. Hence the external 
gravitational potential is given by 

4/ext(f,x)=^^^-^, (10) 

where M, and x,(f) are the mass and positions of the corre¬ 
sponding objects. Regarding the self-gravity of the flow, solv¬ 
ing equation m within is equivalent to computing the 
integral 

Once 4*ext and 4'int are obtained, we can then use equation 0 
to obtain the total gravitational field and subsequently use it 
in both the Navier-Stokes equation and in integrating the tra¬ 
jectories x,(f) for the external objects. Although for our test 
problems we assume that the central object does not move and 
solve the hydrodynamic equations in a fixed reference frame, 
it is trivial to generalize our algorithm to co-moving coordi¬ 
nates. 

There are two different approaches to reducing the above 
three-dimensional formalism to two dimensions, depending 
on the physical problem under study. The first one assumes 
that the density is independent of the vertical coordinate z, i.e., 
Pint(f,r,(^,z) = Pmd(f,T,0). In this case, we define ijjit,r,(j)) = 
4'int(f)T, 0,z) and we are left with a two-dimensional prob¬ 
lem. We obtain the gravitational potential by solving the two- 
dimensional Poisson’s equation 
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(This approach is used in ZEUS-2D, see IS tone & NormanI 
1 19921 for details.) 

The second approach assumes that the vertical structure of 
the density is described by some function Z{r,z) that is inde¬ 
pendent of t and (j), i.e., pi„t{t= S(f,r, ^)Z(r,z). We 
are interested in the gravitational potential on the z = 0 plane, 
so that we solve for ilj{t,r,(j)) = = 0). This is 

a “pseudo-two-dimensional” problem, where the potential is 
not described by a two-dimensional Poisson equation. In or¬ 
der to compute the potential properly, the easiest method, in 
this case, is to integrate directly the equation 


'4’{t,r,4>) = [ ,4l)T.{t,r',(l)')r'dr'd(j)', (13) 

where denotes the two-dimensional computational do¬ 
main and 


,(j)') = -G j 


zy,z')dz' 


y P-+y - ird cosy - (j)')+y 


is the “modified Green’s function” for our problem. 


(14) 


3. NUMERICAL METHODS 

In this section, we describe the numerical methods we use 
to solve th e two classes of s elf-gravity problems. As in the 
first paper llChan et al.ll2005l) . we use pseudo-spectral meth¬ 
ods, in which we expand all functions in series by choosing 
the truncated functions to agree with the approximated func¬ 
tions exactly at the grid points. 

We use the modified Chebyshev collocation method along 
the radial direction and the Fourier collocation method along 
the azimuthal direction. Let N+l and M be the number of 
collocation points in the radial and azimuthal direction, re¬ 
spectively. For every function /(r, (j>), we use the subscript m 
to indicate the Fourier coefficients 

M/2-1 

firy)= E (15) 

m=-Mj2 


for m = 0,1,... ,M/2, where we have set / = 47rGpind. Note 
that, in discrete Fourier transforms, /o and ^re real but 
the other coefficients are complex, so there are, in total, M 
independent equations. We use A,„ to denote the differential 
operator for each m. In the case of the r-direction, the Cheby¬ 
shev polynomials are not eigenfunctions of A„,. Nevertheless, 
because A,,, is time independent and linear, we can write Am 
in matrix representation and pre-compute its inverse A“/ with 
an solver for each m. 

Let ? be the standardized coordinate of Chebyshev polyno- 
mials and r — s[ r) be the mapped (physical) coordinate (see 
IChan et ^120051 for more details). It is well known that the 
spectral Chebyshev derivative in the standardized coordinate 
can be written in matrix representation as 


2^2+1 

6 

... ... 

l-Xj 

li-m \ 

1 (-1)' 

2 \-Xi 

i-iy 

Xi-Xj 

-^i 

2(1-^) 

X-Xj 

1 (-1)'^+'' 

2 1+A,- 

-i(-l)A 

9 (-1)""^ 

Uxj 

2^2+1 

6 / 


(18) 


which is a (A -I- 1) x (A -I- 1) matrix (discussion of 
computing this matrix accur ately can be found in 
iBaltensperger (Q’rummerl 120031) . The derivative in the 
mapped coordinate can be obtained by the chain rule 


d _ dr d 1 d 
dr dr dr dg/drdr' 


(19) 


so that the derivative matrix in the mapped coordinate is given 
by the product 


{dg/dr)\r,^‘^' 


( 20 ) 


With this we are able to write 


and the subscript n to indicate the Chebyshev-Fourier coeffi¬ 
cients 

N M/2-1 

/('•,</>) = £ L fnmTn(r)e^"'\ (16) 

n=Qm=-Ml2 

where f € f-1,1 ] denotes the standardized coordinate (see 
IChan et d for notation). 

We implement a standard two-dimensional ps eudo -spectral 
Poisson solver, which we describe briefly in ^3.11 (a good 
introduction is available in lT^fethenll2000l) . For the direct 
integrator, we describe in ^.3.21 how to gen eralize the com¬ 
pact c ylindrical Green’s function expansion (tCohl & Tohlinei 
1 1999ft with pseudo-spectral methods. Both methods use pre¬ 
computed matrices to speed up the algorithms. The run-time 
computational cost for both methods are of order 0{N^M + 
A/WlogjM). 


3.1. Two-Dimensional Pseudo-Spectral Poisson Solver 

We describe here a two-dimensional, pseudo-spectral Pois¬ 
son solver. We use the fact that the basis polynomials in 
the (p direction, are also eigenfunctions of the operator 
d^/dy in order to split the Poisson’s equation to 


Am'l/’m — 


/ d^ d / 

rdr r^ J ™ 


fm 


(17) 


A, 


m,ij 


^DikDkj 

k 


^ik ^ij 2 

J.2'” ■ 


( 21 ) 


The matrix Am,,; is, of course, singular because a unique so¬ 
lution to tjj,„ does not exist until the boundary conditions are 
given. 

An important case is the vanishing boundary condition 

l/’m(l'min) “ lAm(l'max) — O- Am,,; be the (A— 1) X (A— 1) 
sub-matrix of Am,,/ generated by removing the first and last 
rows and columns. Imposing the vanishing boundary condi¬ 
tion to A“\.^ is then equivalent to finding the inverse of the 

Am,,; and filling back the first and last rows and columns with 
zeros, i.e.. 


/o 

0 • 

• 0 


0 

A-'.. 

0 

0 

rn,ij 

0 

\0 

0 • 

• 0 

0/ 


We can then calculate 

4?Vo = E 

j=0 


( 22 ) 


(23) 
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where the superscript (0) indicates that the solutions sat¬ 
isfy the homogeneous boundary conditions for all m. Tak¬ 
ing the inverse Fourier transform along the (/(-direction, we 
then obtain the solution which satisfies the vanish¬ 

ing boundary conditions. 

In order to apply more generic boundary conditions, we first 

look for solutions ipm\r) that satisfy the homogeneous differ¬ 
ential equation 

^ 0 (24) 

with the proper boundary conditions. Then it is clear that the 
sum +ipm^ satisfies 

+ =/m (25) 


with the same boundary conditions. The solutions to equa¬ 
tion are 


j , m = 0 

^ + , m= 1 , 2 ,..., M/2 


(26) 


where Cm and D„, are complex constants, which are fixed by 
the boundary conditions. In general, these boundary condi¬ 
tions depend on azimuth, i.e., 

=An(</>), (27) 

'/’('■max,/') =/3out (</>)• (28) 

We first take the Fourier transform of /3in(</') and /3out(</'), we 
solve for and Dm using 


Pin,m — j 

f Cfn In Y + DfYi !) ^ — 0 

+ /2,. 

..,M/2 

(29) 

/5out,m — j 

f Cm ^min , HI = 0 

lC’mCn + ^"F-;;Tn,'M= 1,2,. 

..,M/2, 

(30) 

we add these terms to /m ^, and finally take the inverse Eourier 


transform. 

We summarize the two-dimensional pseudo-spectral Pois¬ 
son solver using the following steps: 


1. We take the Fourier transform of each physical quantity 
/ along the (/-direction and obtain fm by a fast Fourier 
transform, which is of order ff{NM\og 2 M). 

2. We then compute the matrix product J23> for each m, 
which is of order 


3. We compute /3m by a fast Fourier transform, which is of 
order i^(Mlog 2 M). 

4. We then solve for Cm and Dm using equations (I29> and 
(EOj, which is of order 

5. We impose the boundary conditions using equa¬ 
tion for each m, and obtain 'tpm{rk), which is of 
order 

6. Finally, we take the inverse transform of ijjm to obtain 
the potential /;, which is of order i^(A'Mlog 2 M). 

Therefore, the overall computational cost is + 

NM\og2M). 


3.2. Two-Dimensional Gravity Integrator 

As described in ^ we assume in this case that the vertical 
structure of the density is time-independent and has cylindri¬ 
cal symmetry. Hence we can write 

Pint(f,r,(/),z) = S(f,r,(/))Z(r,z), (31) 

where we normalize Z{r,z) along the z direction to unity, i.e., 

J Z{r,z)dz=l, (32) 

so that S has the physical meaning of a height-integrated den¬ 
sity. 

We present here a new numerical method to compute the 
integral 113) efficiently and acc urately. Our method is moti¬ 
vated bv lcohl & Tohlin j 119991) . who showed that the three- 
dimensional Green’s function in cylindrical coordinates can 
be expanded in the compact form 






Sm-l/2(x), 


(33) 


where x = [r^-t-r''^-t-{z—z!)'^]/2rr' and 2 m-i/ 2 (x) the half¬ 
integer degree Legendre functions. Note that the </ and (/' de¬ 
pendence appears only in the exponential ). This can 

be done because the Green’s function itself can be written as 
afunction of A (/ = (/-(/' so that 2m-i/2(x)/''’")/'^ simply 
the Fourier transform of |x-x'|“^ with respect to A(/). This 
property is still true for the modified Green’s function because 
Z is /-independent: 

^#(r, (/;/,(/')= ^#(r,r',A(/)). (34) 

Expanding it in Fourier series, i.e.. 


^(r,(/;/,(/.')= £ (35) 

m=—oo 

and substituting it into equation O, we obtain 

p oo 

tp{r,(j))= / Y, ^m{r,r')e‘"'^^~'^'y{r',(j)'ydr'd(j}'. (36) 

' m=-oo 

Note that the integration over (/', 

/""uf(/.'S(r',(/')e-™‘^' = 2TTEm{r'), (37) 

Jo 

is the Fourier transform of E(r',(/)') in the azimuthal direction. 
Therefore, we can rewrite equation as 

i’{r,4>)= Y / %i{r,r')2TrY.m{r')e‘"^^r'dr'. (38) 

m=—<x> ^min 

Taking the Fourier transform of the whole equation again with 
respect to /, we obtain a one-dimensional integral for each 
value of m 

^ /■7'max 

= / ^m{r,r')2TrEm{r')r'dr'. (39) 

^min 

We now expand Sm('‘0 in Chebyshev polynomials, i.e.. 


tmy) = Y'^n,nUf). (40) 

n=0 

Recalling that f' is the standardized coordinate, we let F = 
g(F) be some coordinate mapping. Equation 139) can then be 
rewritten as 

^ ffmax 

^ / %t{r,r')2TTTn{r')r'dr'. (41) 

H=0 
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Fig. 1 .— Plots of the quantity J^i„{r';r) for the first six values of m, when N = 33, M = 64, and Z{r,z) = S{z). 


Because the orthogonality condition of Chebyshev polynomi¬ 
als has a weight function (1 i.e., 


{Ti 




w)un tt ^ 

—dr = — (5/„, 

yj _ ^/2 Cn 


(42) 


where cq = 2 and c„ = 1 for n = 1,2,..., it is not difficult to 
see that, by defining 


ym(r;r) =^m [r,g{r)] Wl-r’^ (43) 


we can write 


J. r..\ 

Wtnv) — ^ ^nm I r - -p^ dr 

„=0 J-1 Vl-f'2 


n=0 

,1=0 


(44) 


where J^nm{r) is the «-th Chebyshev coefficient of J^m{i^',r) 
with respect to r'. 

Up to this point, although the function Y,{r',(jJ) has been 
written in Chebyshev-Fourier series, we have assumed that 
the series is infinite and hence the formalism is exact. In order 
to perform it numerically, we have to discretize equation il. 

For the azimuthal direction, which is periodic, we simply 
replace the continuous Fourier transform by a discrete Fourier 
transform and use the proper normalization, which does not 
change equation (EU but only limits the index m to the range 
-M12 to M/2 -1. For the r-direction, we truncate the infinite 
sum at a finite number of terms, N. Therefore, equation ii 
becomes 

Wm{rk) = V -, (45) 

n=0 ‘ 2 « 

where denotes the collocation points in the r-direction. 

Once we obtain J^nmifk) for n = Q,l,...,N, the calcula¬ 
tion of the gravitational potential is trivial. Flowever, note that 
2m-i/2(x) is related to the complete elliptic integral, which is 


singular when r' = r^-. In general, this difficulty arises because 
of the singularity of |x-x' |“*. In order to avoid the singularity, 
we use the “Chebyshev-roots grid”: 


fj — cos 


(2/+l)7r 

2N 


0<j<N 


(46) 


so that J^m(dj',rk) is finite for all m and j. In Figure[2 we plot 

; r) for the first six values of m for N — 33, M = 64, and 
Z{r,z) = (5(z), which is the case of a very thin disk. 

The Chebyshev coefficients J^nm are now well defined as 


N-l _ N-l 

j) — ^ ‘^nmTni/j) — ^ 
n=0 n=0 


Tm{2j+ 1) 


2N 


(47) 


Note that J^nm cannot be computed by a standard discrete co¬ 
sine transform. It is related to the discrete Fourier transform 
with different parity properties. This non-standard cosine 
tra nsform is known as type -II discrete cosine transform (see, 
e.g. lFrigo & .lohnson l2003i p.41). The computational order is 
still ff{Nlog 2 N). The only potential inconsistency about this 
method is that there are only N collocation points, instead of 
N+1. Nevertheless, because the Chebyshev coefficients con¬ 
verge exponentially for smooth functions, Eai„, is assumed to 
be very small. The artificial vanishing of is, therefore, 
negligible in the final solution. 

In our implementation, we pre-compute the operator 
J^km{fk) to speed up the algorithm. The steps to setup the 
gravity integrator can then be summarized by the following: 


1. We first take the discrete Fourier transform of E(r^, (/;■) 
along the (/-direction and obtain Yl„{ri), which is of 
order ^(NMlogjM). 


2. We then calculate the matrix product iprnirk) = 
Yli=o'^m{ri)./mi[rk)/Cn for each m, which is of order 
0{N^M). 


3. Finally, we take the inverse Fourier transform of '4>m{rk) 
and obtain the restricted potential '0(r^, c/j), which is of 
order (^(NMlogjM). 

The overall computational cost is ff{N'^M+NM\og 2 M)- 
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3.3. Coupling to Hydrodynamic Equations 

We evo lve _ die hydrodynamic equations following 

IChan et alJ ( [200l . i.e., using a low-storage, third-order 
explicit Runge-Kutta method. The only difference is that, at 
the beginning of every time step, we update the gravitational 
potential. 


4. CODE VERIFICATION 

In IChan et al ] j200l , we verified the hydrodynamic part of 
our algorithm. In this second paper, we will first carry out a 
few time-independent tests for both the Poisson solver and 
the direct gravity integrator. We use for each test an ana¬ 
lytic density-potential pair and show that the numerical so¬ 
lutions agree with the analytic expressions. We also perform 
some time-dependent tests, which demonstrate that the grav¬ 
ity solver couples to the hydrodynamic equation correctly. 


4.1. Time Independent Tests of Poisson Solver 

The simplest way to test the Poisson solver is to com¬ 
pare some numerical potential i/^num obtained by the Pois¬ 
son solver to its analytical expression '0ana. The tests in 
this subsection are done in the computational domain = 
[0.2,1.8] X [-tTjTt). We consider a very simple potential 




2 1,0. 0.0648 \ 

r^-a\ 1.82r-) 


sinr 


(48) 


where a is some arbitrary parameter. We chose this potential 
because its Laplacian takes the form 

/ana(r, f) = V^'f/'ana = sin (49) 

which is independent of a. Hence, the parameter a forces our 
solution to satisfy different boundary conditions. 

We choose three different boundary conditions to test the 
Poisson solver. The first one is 

I l/'ana(rmin,<^) = jr^jnSinc^ ^ 

which corresponds to cr = 0. For the second one, we choose 
the vanishing boundary condition 

f V^ana(Tmin5 — 0 (51) 

which corresponds to cr = 1. Finally, we also choose 

/ V^ana(rmin; — 3 (I'niin ~ ^*^^^rnin4" 0- 1296/rj.nij.i) sin (ji 
\ V^ana(l'max; ^ ~ 3.64rniax3" 0.1296/riTiax) sin 

(52) 

which corresponds to cr = 2. In Figure|2] we show the contour 
of /ana, V^num, and the difference, e = ji/’num-'/anaj, for each 
choice of cr. The contour lines for finum show clearly how 
the boundary conditions change the solution. The numerical 
results agree very well with the analytical solution and the dif¬ 
ference is only of order 10“*^, which is the machine accuracy. 


4.2. Tme Independent Tests of Gravity Integrator 

To test the gravity integrator, we consider an infinitesimally 
thin disk, with an exponentially decaying surface density, i.e., 

= s:{r,(j))5(z) = (53) 

The corresponding potential on the z = 0 plane is given by 

=-7iGEor[/o(y)/:i(y)-/i(y)/:o(y)] (54) 


where y = r/2cr (iBinnev & Tremain3ll987h . The functions 
In{y) and Kn{y) are the modified Bessel functions of the first 
and second kinds, respectively. The modified Green’s func¬ 
tion for this problem is simply 

-G 


^{r,(j),r',(j)') = 


(55) 


-t r'2 _ 2rr' cos (/ - /') 

In order to test our gravity integrator for non-axisymmetric 
problems, we place three disks in the computational domain 
with each disk centered at some (r,-,/,). Let 


Ri = r^ + rf — 2rriCOs{(f>— (fj) 


(56) 


be the distance between {ri,cj)i) and (r,/). Normalizing the 
total mass for each disk to unity, we have 

1 / R, 




){rA) = 


2'Ka'^ 


exp 


(57) 


The corresponding potential is 


i’{r,,,Pi){T4')=--yi[Hyi)K\{yi)-I\{yi)KQ{yi)\ , (58) 

where y,- = R,/2 ct. 

Because the potential is singular, we put the center of 
each disk off grid. Specifically, we choose the coordinates 
(0.9,7r/4), (0.9 ,tt), and (l,-7r/3). We also change the total 
mass of each disk, in order to avoid cancellation of errors by 
symmetry, by multiplying them by different constants so that 

Eana(r,/) = E(o.9,.n./4) (r,/)-t -£( 0 , 9 ^,^.) (r,/) -F2E(i__.,r/3) (f/) ■ 

(59) 

We repeated this test for different choices of a, namely, 
cr = 0.025, cr = 0.05, and cr = 0.1. In Figure 0 we show 
the grayscale plots for Sana, '/num, and the fractional error 
e = |'/num/'0ana“ 1|, for each value of cr. The contour lines 
in the error plots show that, at the smooth density region, the 
error is of order 10“^. The maximum errors appear around the 
density peak, where the analytic potential is singular. 

Besides being able to handle a (5-function in the z-direction, 
our gravity integrator is able to solve problems with other 
time-independent and (()-independent vertical structures. For 
example, we consider the Gaussian sphere given by 

= tzTTTTTTT exp 




,2+z2l 


2(7^ 


(60) 


which is normalized so that the total mass is M. Here ct is a 
parameter that controls the spread of the mass in each sphere. 
As cr ^ 0, the Gaussian sphere approaches a point source. 
If we have a collection of such spheres with the same a, the 
vertical structure can be factored out of the equations. 

Using Gauss’ law, we find that the potential on the z = 0 
plane is given by 


ip{r, (p) = —erf 
r 

where erf(x) is the error function 

2 


(V2a) 


erf(x) = —= [ e ^^df. 


(61) 


(62) 


For a collection of Gaussian spheres centered on the z = 0 
plane with the same a, we can factor out the z-dependence, 
i.e. 

p{r,<P,z) = Y,P{ri,^,)iT(t>,z) =Z{r,z)Y,^{ri,4,i){r,<P) , (63) 
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Fig. 2.— Grayscale plots of the analyt ical f unction /ana, the numerical potential and the difference e between the numerical solution and the analytical 
solution for the test problem discussed in EH In all plots, darker shades correspond to larger magnitude. 


where we have used the same notation as in the previous sub¬ 
section to indicate that the center of each sphere is located at 
(r, , (/),). The normalized vertical structure is given by 

The surface density for each Gaussian sphere is, therefore, 

1 f \ 

S(.A)(T<^) = ^exp(^-^j, (65) 

where Rj is defined in equation (|56}. 

The modified Green’s function for our problem is 


r,(l)-y ,(!)')= j 


-Z{r',z')dz' 


-<» yr^ + r'^ — 2rr'cos((j)— (/)')+ z'^ 
e‘^'^^Ko{R^/4) 


V2' 


( 66 ) 


TTCr^ 


where Kq denotes the modified Bessel function of the second 
kind and we have set G = 1 for simplicity. We use the modi¬ 
fied Green’s function described in equation ( l66l to obtain the 


correspond ing Then we follow the algorithm de¬ 

scribed in 3.3.2l to perform the calculation for the gravitational 
field for three Gaussian spheres located at (r, = (1,0), 

(0.75,27r/3), and (1.25,-27r/3) with a total mass of 2, 1/2, 
and 1, respectively, i.e., 

Sana(6 />) = 2S(l,o) (r, (j)) + iS{o.9,37r/4) />) + ^{l-7v/2) </) 

(67) 

In Figure 13 we show the numerical solutions for a = 0.05, 
tr = 0.1, and a = 0.2. The maximum numerical error is less 
than 0.5%. 

4.3. Free Fall of a Dust Ring Under Self-Gravity 

In order to perform a dynamic test of our spectral self¬ 
gravity solver, we follow the free-fall of a self-gravitating dust 
ring. The initial conditio n of this problem is the same as in 
§4.1 of iChan et alJ Jliol . i.e., 

P{ry,z) =pir,<p) =exp[-20(r-l)^], (68) 

and the initial velocity is zero. The only difference is that, 
this time, we explicitly calculate the self-gravity of the flow 
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Fig. 3.— Grayscale plots of the density S, the potent ial -0 , and the difference e between the analytical and the numerical values of the gravitational potential 
for a selection of infinitesimally thin exponential disks (' ^4.2i . The top panel corresponds to cr = 0.025, the middle panel corresponds to (T = 0.05, and the bottom 
panel corresponds to cr = 0.1. 


instead of using the gravitational field provided by a central 
object. We set the gravitational constant to G = 1 and let the 
ring free fall. We neglect pressure and viscosity and set the 
resolution to 513 x 64 in order to resolve the sharp peak in the 
density prohle that appears at late times. 

In order to test our implementation of the spectral method 
and in the absence of an analytic solution to the problem we 
wrote a very simple N-body algorithm. We placed 100,000 
particles in the computational domain based on the initial den¬ 
sity of the ring. Then we computed the gravitational interac¬ 
tions between each particle pair directly and integrated their 
trajectories. 

In Figure|3 we plot the density prohles of the numerical so¬ 
lution of the free falling dust ring at (/) = 0 using the spectral 
algorithm and the N-body method. Different lines represent 
the solution from our pseudo-spectral algorithm at different 
times and the open circles represent the solution from the par¬ 
ticle method. The two solutions agree very well. 

4.4. Orbiting Cylinders 


We adopt this test from iFrver et alJ ll2005h . who simulated 
two spheres orbiting around each other with smooth particle 
hydrodynamics. Because our algorithm is two-dimensional, 
we assume that all variable are independent of z and mod¬ 
ify the problem to two inhnite cylinders orbiting around each 
other. 

Because the region outside the cylinders has very low den¬ 
sity, small fluctuations of the solution in this low-density re¬ 
gion result in a large error in the pressure force. In order to 
resolve the numerical instability that arises from this effect, 
we evolve the whole set of hydrodynamic equations including 
the energy equation and introduce a small background term 
in both the continuity equation and the energy equation. This 
background needs to be small compared to the physical prop¬ 
erties of the cylinders so that it does not affect the orbital mo¬ 
tion signihcantly, while at the same time it needs to be large 
enough so it can screen out the errors at the low density region 
in our algorithm. Then, a strong spectral hlter is applied to the 
velocity field to reduce numerical instabilities. The spectral 
biters for density and energy, on the other hand, are relatively 
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Fig. 4.— Grayscale plots of the density S, the poten tial ij ), and the difference e between the analytical and the numerical value of the gravitational potential 
of a collection of three-dimensional Gaussian spheres ( ^4.21 . The top panel corresponds to cr = 0.05, the middle panel corresponds to cr = 0.1, and the bottom 
panel corresponds to cj = 0.2. 



Fig. 5.— The radial density profile of the numerical solution for a free- 
falling dust ring at ^ = 0. The solid lines represent the solution using our 
spectral algorithm and the open circles represent the solution using a simple 
N-body code. 


We consider a Gaussian density distribution for each cylin¬ 
der 


E(r,0) = 


27r(T^ 


(69) 


In order to find the initial condition for the energy equation, 
we first need to solve for the gravitational potential using the 
Neumann boundary condition dtp/dr — 0 at the origin. The 
second boundary condition is not important here because it 
only shifts the potential by a constant, which does not affect 
the gravitational field. Setting the arbitrary constant for the 
second boundary condition to zero, we obtain 


tP = G 


21nr-Ei - 


2cr^ 


(70) 


where Ei(.v) is the exponential integral function defined by 


Ei(x) = - 


f 


-dx' 


(71) 


The corresponding gravitational field of this potential is 
2G 




weak in order to conserve mass and total energy. 


r 


(72) 
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Fig. 6.— Grayscale plots of the density E of the two orbiting cylinders and vectors plots of the velocity fields at different times in the simulation. In all plots, 
darker colors correspond to larger magnitudes. 



Fig. 7.— The total angular momentum of the orbiting cylinders (normal¬ 
ized initially to unity) as function of time. The angular momentum is con¬ 
served to better than 2% over t — 100, which corresponds to about 16 periods 
or ~ 5 X 10*^ timesteps. 

Therefore, we require the pressure to be 

P{rA) = j ^8 /la^)] 

(73) 

in order to balance self-gravity. Note that we again set the 
integration constant to zero in the pressure. The modified 
Green’s function for this problem is 

^^(r,= Gin [r^-l-/^-2r/cos((^-(/)^)] . (74) 

In our simulation we have two Gaussian cylinders centered 
at (r,■,(/),■). Therefore, we set 

S(r,*)(--.« = - 5 -^ (75) 


and 

^ [Ei{-Rj/a^)-Ei{-Rj/2a^)] , (76) 

where /?, is defined in equation ( I5SI . Because the pressure 
is singular when Ri = 0, we choose the initial position of the 
cylinders {ri,(f)i) off grid. In particular, we also choose the 
background density such that the area between the cylinders 
contains only 1% of the total mass in the whole system. Be¬ 
cause the area of our domain is 7r(1.8^-0.2^) = 3.27r, we set 


E{r,(j)) =0.99 E(i 10-3)(r,(/))-!-E(i ,,.,10-3)(r,(/() 


and 


7^(6 «<>)= 7 ^( 1 , 10 - 3 ) (6 </()+^’( 1 ,^+ 10 - 3 ) 


0.02 

3^' 


0.02 

'T2^ 

(77) 


(78) 


The initial energy is then chosen using the equation E = 3P/2 
that is based on an ideal gas law. Setting G = 1, the angular 
velocity of the cylinder needed to balance the gravitational 
force is equal to unity. Hence, we set the initial velocities to 


Vr = 0, = r. (79) 

We use cr = 0.1 and evolve the simulation to a dimension¬ 
less time t = 100, which is equal to about 16 complete orbital 
rotations. In Figure we show the gray-scale plots of the 
density as well as the velocity field for different times in the 
simulation. In Figure 0 we plot the angular momentum as 
a function of time. Our algorithm is able to conserve angu¬ 
lar momentum to better than 2% for up to about 16 complete 
orbiting motion, which corresponds to ^ 5 x 10^ timesteps. 
The 2% dissipation of angular momentum is mostly due to 
the strong spectral filtering, which is equivalent to artificial 
viscosity, and the deformation of the two cylinders during the 
simulation. 
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Fig. 8. — The minimum value of the Toomre parameter min(2) anywhere 
in the disk, for different values of the polytropic constant K. The three dif¬ 
ferent lines correspond to different values of the polytropic index F. When 
min{Q) < 1, i.e. in the shaded region, the linear analysis suggests that the 
disk is unstable. Filled circles denote simulations in which the disk was un¬ 
stable, whereas open circles denote simulations in which the disk was stable. 

5. THE STABILITY OF SELF-GRAVITATING HYDRODYNAMIC 
DISKS 


Pseudo-spectral methods effectively evolve in time the hy¬ 
drodynamic equations for each mode independently. For this 
reason, they are ideal tools for confirming the results of lin¬ 
ear mode analysis and for extending them to the non-linear 
regime. As a case study, we address here numerically the sta¬ 
bility of self-gravitating infinitesimally-thin disks. 

We consider self-gravitating gas disks with a polytropic 
equation of state 

P = KJ:^, (80) 


where K is the polytropic constant and F is the poly¬ 
tropic index. The disks ar e locally stable to all 
non-axisymmetric perturbations jGoldreich & Lvnden-Belll 
Il96-5aibt Llulian & Toomrell1966h . For axisymmetric pertur¬ 
bations, linear-mode analysis of the hydrodynamic equatio ns 
results in the dispersion relation llBinney & Tremaine! 1987t) 

UJ^ = K^-2TTG^k\+k^cl, (81) 


where we use Cs to denote the sound speed of the gas, k to 
denote the epicyclic frequency in the disk, and uj and k to de¬ 
note the mode frequency and wayenumber, respectiyely. The 
disk is locally stable if 0. The condition for all modes 
^tt^b^st^e i^ive^by Toomre’s stability criterion (ISafronovI 


Q = 


CsK 

ttGS 


> 1 . 


(82) 


Although criterion (m is only a local condition, it is also 
sufficient for global stability of axisymmetric modes with 
kr ^ 1. At the critical yalue 2=1, the disk can be glob¬ 
ally unstable to non-axisymmetric modes. As Q increases, the 
disk becomes globally stable to all non-axisymmetric modes. 
Howeyer, for yery large yalues of Q, pressure dominates the 
self-grayity effects and the disk becomes globally unstable to 
non-axisymmetric perturbations again. Our goal in this sec¬ 
tion is to study numerically the stability of self-grayitating 
disks in situations where the unstable wayelengths are larger 
than the characteristic length scales in the disks. 

We simulate the time e volution of disks that obey Plum¬ 
mer’s density model (see iRinnev & Tremair^ l1987l chap¬ 
ter 2) 


E(r) 


Mja 

2TT{r'^ + a^)^^^ ’ 


(83) 


where stands for the initial disk mass and cr is a parameter 
that controls the central concentration of the disk density. The 
corresponding gravitational potential (at the z = 0 plane) is 
given by 


V’(t) 


GMd 


(84) 


We set the initial radial velocities to zero, i.e., Vr = 0, and 
choose the azimuthal velocities, v^, so as to balance the pres¬ 
sure and self-gravity of the flow, i.e.. 


_ 27rGEr 3KrT,^-'r 
r a r^ + a^ 


(85) 


Although Plummer’s model assumes that the density is a 6- 
function along the z-direction, its corresponding scale height 
H = rcs /is given by 


// _ Cs _ / 27rGr2 _ 3r^ \ 

r Vcj) 


( 86 ) 


We choose 129 x 64 grid points and perform the simulations 
in the domain [0.005,5] x [-7r,7r) with cr = 1. We set the total 
disk mass to = 1. 

Figure 13 a parameter study of the dependence of the stabil¬ 
ity of the disk on the polytropic index F and the constant K. 
We plot the minimum value of the Toomre parameter min(2) 
anywhere in the disk against K. The three different lines cor¬ 
respond to different values of the polytropic index F. Accord¬ 
ing to linear analysis, when min(2) < 1, there is some region 
in the disk where the disk becomes unstable. We present the 
results of our simulation on the same plot, with filled circles 
denoting the solutions that become unstable and with open 
circles the solutions that remained stable during the course of 
the simulations. 

For small values of the polytropic index (F < 4/3) or for 
small values of the polytropic constant {K < 0.2), our numer¬ 
ical simulations agree with the predictions of the linear mode 
analysis. The small disagreements near the Q = 1 separatrix 
come from two facts: that, (i) we have a finite domain of inte¬ 
gration on which we have to impose boundary conditions, and 
(ii) we cannot calculate the evolution of the instabilities and 
the perturbations in the gravitational held from matter outside 
the boundaries. Regarding point (i), we have used the ab¬ 
sorbing bo undary condition discussed in the first paper in the 
series Isee lChan^ al.l2005l) . which reduce the reflection sig¬ 
nificantly. However, a very small but finite amount of the en¬ 
ergy of the wave is still trapped in the domain, which causes 
the instability. As a result, when the gravitational potential 
deviates from Plummer’s model because of the finite size of 
the domain, the disk close to the outer boundaries becomes 
unstable. 

A striking difference, however, appears when the polytropic 
index and the polytropic constant become larger: our numeri¬ 
cal solutions are always stable, contrary to the predictions of 
the linear mode analysis. The reason for this disagreement be¬ 
tween the linear mode analysis and our numerical simulations 
becomes evident in Figure|3 where we plot Toomre’s param¬ 
eter as function of radius r for different values of the poly¬ 
tropic constant K. When < 0, there is no initial value of the 
azimuthal velocity that satisfies the steady-state equation 
for radial force balance; we shade this region using dark gray. 
For 0 < 2^ < the disk is unstable according to Toomre’s 
criterion; we shade this region using light gray. Finally, when 
2 "> 1 the disk is linearly stable; we leave this region white. 
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Fig. 9.— The Toomre parameter as a function of radius r for solutions with different values of the polytropic constant K. The dark gray region corresponds 
to < 0, for which no steady-state rotating solution exists; the light gray region corresponds to 0 < Q- < I which is Toomre unstable; the white region 
corresponds to > 1, which is Toomre stable. The contour lines correspond to the size of the minimum unstable wavelength according to the linear mode 
analysis. 


We also plot in Figurecontours that correspond to the min¬ 
imum unstable wavelength, in the unstable region, computed 
by 

^=W = —(l + vT^) (87) 

^min ' 

Figurel^shows that when F > 4/3 and for large values of 
the polytropic constant K, the minimum unstable wavelength 
is larger than the extent of the region in which the disk is un¬ 
stable. For example, for F = 4.5/3 and for K G [0.42,0.56], 
the unstable region only extends to r G [0,1.5]. However, 
in that same region, the minimum unstable wavelength is 
Amin ~ 3, in the appropriate units. Because the unstable wave¬ 
lengths are larger than the characteristic length scale in the 
system (which can be defined here as the extent of the region 
in which 0 < Q < 1), the approximations involved in the linear 
mode analysis are not justified. Indeed, the numerical simula¬ 
tions show that the unstable modes cannot grow and the disk 
is stable. 

6. CONCLUSIONS 

Problems that involve self-gravity are usually time- 
consuming tasks in computational physics. In standard fi¬ 
nite difference methods, Poisson’s equation is solved as the 
steady state of a diffusion equation (using relaxation and 
over-relaxation methods) and has to be recalculated together 
with the hydrodynamic equations at every timestep. Al¬ 
though hybrid algorithms have been developed which use 
high-order methods to solve Poisson’s equation, finite differ¬ 
ence schemes are still used for evolving the hydrodynamic 


equations. The resulting inconsistency in the order of differ¬ 
encing the hydrodynamic equations and Poisson’s equations 
either significantly reduces the accuracy of the Poisson solver, 
or requires the extra complication of interpolating the grav¬ 
itational filed at each grid point. More importantly, existing 
two-dimensional hybrid algorithms can only address a limited 
number of self-gravitating problems because the solutions to 
the two-dimensional and three-dimensional Poisson’s equa¬ 
tion are fundamentally different (see ffland Cl¬ 
in this paper, we presented two different approaches to us¬ 
ing pseudo-spectral methods to solve self-gravity problems. 
For the first approach, we described the implementation of 
a standard pseudo-spectral Poisson solver that solves the two- 
dimensional Poisson’s equations to machine accuracy. Instead 
of solving a diffusion equation, spectral methods allowed us 
to invert the Poisson operator in spectral space, making the 
algorithm fast and accurate. For the second approach, we in¬ 
vestigated a fast gravity integrator for disks-like flows with 
known, time-independent vertical structures. This algorithm 
allowed us to study the evolution of flows with finite, but 
not only infinitesimal, thickness. This improvement allows 
two-dimensional algorithms to solve a whole new class of 
problems. We demonstrated here the ability of our algorithm 
to compute properly and efficiently the gravitational poten¬ 
tial of flattened flows using different test problems. Even for 
the high resolution that corresponds to 129 x 256 collocation 
points, our Poisson solvers and integrator use less than 10% 
of the total computational time. 

We also explored how to extend Toomre’s stability criterion 
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to self-gravitating disks for which the characteristic properties 
of the flows change over a length scale that is shorter than 
the minimum unstable wavelength. Based on our simulations, 
we And that for Plummer’s density model, if the disk is hot 
and has a polytropic index T > 4/3, all oscillatory modes in 
the disk are stable, contrary to the predictions of the analytic 
calculation. 


We thank an anonymous referee for constructive comments 
that improved the clarity of our paper. C.-K. C. and D. P. 
also acknowledge support from the NASA ATP grant NAG- 
513374. 
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